Boolean Network Approach to Negative Feedback 
Loops of the p53 Pathways: Synchronized 
Dynamics and Stochastic Limit Cycles 

Hao Ge* Min Qian t 

April 15, 2009 



Abstract 

Deterministic and stochastic Boolean network models are build for the dy- 
namics of negative feedback loops of the p53 pathways. It is shown that the 
main function of the negative feedback in the p53 pathways is to keep p53 at 
a low steady state level, and each sequence of protein states in the negative 
feedback loops, is globally attracted to a closed cycle of the p53 dynamics 
after being perturbed by outside signal (e.g. DNA damage). Our theoreti- 
cal and numerical studies show that both the biological stationary state and 
the biological oscillation after being perturbed are stable for a wide range of 
noise level. Applying the mathematical circulation theory of Markov chains, 
we investigate their stochastic synchronized dynamics and by comparing the 
network dynamics of the stochastic model with its corresponding determin- 
istic network counterpart, a dominant circulation in the stochastic model is 
the natural generalization of the deterministic limit cycle in the deterministic 
system. Moreover, the period of the main peak in the power spectrum, which 
is in common use to characterize the synchronized dynamics, perfectly corre- 
sponds to the number of states in the main cycle with dominant circulation. 
Such a large separation in the magnitude of the circulations, between a dom- 
inant, main cycle and the rest, gives rise to the stochastic synchronization 
phenomenon. 



* Centre for Computational System Biology and Mathematical Department, Fudan University, 

Shanghai 200433, P.R.China; email: edmund_ge@tom.com 

'''School of Mathematical Sciences, Peking University, Beijing 100871, P.R.China 



1 



KEY WORDS: Hopfield network; Stochastic Boolean network; Boltzmann 
machine; circulation; stochastic limit cycles; synchronization; p53 protein; 

I Introduction 

Thousands of papers have been reported in the field of p53 protein during the last 
three decades since 1979, including both the experimental and theoretical analysis. 
p53, the tumor suppressor, transcriptionally activates Mdm2, which in turn targets 
p53 for degradation (Piette et al. 1997, Prives et al. 1998, Vogelstein et al. 2000, 
Momand et al. 2000, Michael et al. 2003), keeping p53 at a low steady state level. 
The low steady state of p53 is essential for normal cell prohferation because p53 
induces either cell cycle arrest or programmed cell death. The concentration of p53 
increases in response to stress signals, such as DNA damage, inducing a transition 
to oscillations of p53 level. Following stress signal, p53 also activates transcription 
of several hundred genes, which are involved in the program of cell cycle arrest, 
apoptosis, cellular senescence and DNA repair. Therefore, it is important to note 
that many additional proteins interact with p53 and Mdm2, so that a number of 
positive and negative autoregulatory feedback loops acting upon the p53 response 
are embedded inside a network with a few additional interactions (Harris et al. 
2005). 

Several mathematical models have been proposed to explain the damped oscilla- 
tions of p53, either in cell population or in a single-cell, most of which are determin- 
istic model of ordinary differential equations (Mihalas et al. 2000, Lev Bar-Or et al. 
2000, Monk et al. 2003, Ma et al. 2003, Ciliberto et al. 2005). On the other hand, 
the stochastic nature of p53 dynamics has also been observed and caused more and 
more interest nowadays (Ma et al. 2003). 

However, in many cellular biochemical modelling, a detailed model, whether 
the deterministic one based on mass-action law, Michaehs-Menten kinetics and Hill 
function, or the stochastic molecular number-based CME(chemical master equa- 
tion), is not warranted because of a lack of enough quantitative experimental data. 
Alternatively, one can develop discrete state network model, i.e., deterministic and 
stochastic Boolean networks, of a complex biological system using the available in- 
formation on the activation and repression from one signaling molecules to another, 
e.g., the kind of signahng wiring diagram as the Fig. 1 in (Li et al. 2004), where 
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Li et. al. have developed a discrete Boolean model by applying the approach of 
Hopfield for neural networks for yeast cell-cycle regulatory network with 11 nodes. 
The main results of (Li et al. 2004) are that the network is both dynamically and 
structurally stable. The biological steady state, known as the Gl phase of a cell 
cycle, is a global attractor of the dynamics; the biological pathway, i.e., the return- 
ing to Gl phase after perturbation, is a globally attracting dynamic trajectory. The 
deterministic Boolean model has been further extended to incorporate stochastic 
dynamics (Zhang et al.2006). They found that both the biological steady state and 
the biological pathway are well preserved under a wide range of noise level. Fur- 
thermore, recently we investigated the synchronized dynamics and nonequilibrium 
steady states in the stochastic yeast cell-cycle network (Ge et al. 2007), applying 
the mathematical circulation theory of Markov chains( Jiang et al. 2004). 

In the present paper, deterministic and stochastic Boolean network models are 
build for the dynamics of negative feedback loops of the p53 pathways. It is shown 
that the main function of the negative feedback in the p53 pathways is to keep p53 at 
a low steady state level, and each sequence of protein states in the negative feedback 
loops, is globally attracted to a closed cycle of the p53 dynamics after being per- 
turbed by outside signal (e.g. DNA damage). Our theoretical and numerical studies 
show that both the biological stationary state and the biological oscillation after 
being perturbed are stable for a wide range of noise level. Applying the mathemati- 
cal circulation theory of Markov chains, we investigate their stochastic synchronized 
dynamics and by comparing the network dynamics of the stochastic model with 
its corresponding deterministic network counterpart, a dominant circulation in the 
stochastic model is seen to be the natural generalization of the deterministic limit 
cycle in the deterministic system. 

It is shown that a large separation in the magnitude of the circulation, between 
a dominant main cycle and the rest, gives rise to the stochastic synchronization 
phenomenon and the stochastic global attractive behavior, and moreover the power 
spectrum of the trajectory has a main peak, whose period converges just to the num- 
ber of states in the dominant cycle. Furthermore, the net circulation of the dominant 
cycle increases monotonically with the noise-strength parameter /3, approaching its 
deterministic limit; while the circulation of all the other cycles approaches zero very 
fast when P is quite large. Together, these observations provide a clear picture of the 
nature of the synchronization and stochastic limit cycles in a stochastic network in 
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terms of the probabilistic circulation of NESS (nonequilibrium steady states) (Jiang 
et al. 2004). 

For the completeness of the work, in supporting information, we give a theoretical 
sketch of some relevant results on biological networks, including a classification of 
the deterministic and stochastic Boolean networks and their correspondence. Also a 
short introduction of the mathematical theory of stochastic circulation for Markov 
chains is introduced and applied to deterministic and stochastic networks. It is 
shown that the stochastic Boolean network is reversible if and only if the matrix T 
in the model is symmetric, and the net NESS circulation is strictly positive as long 
as the probabihty of the directed cycle is larger than that of its reversal. 



II Boolean network approach 

Since the influential work of J.J. Hopfield in 1980s' (Hopfield 1982, Hopfield 1984), 
the deterministic Boolean (Hopfield) network has been applied to various fields of 
sciences. Amit (Amit 1989) has introduced a temperature-like parameter /? that 
characterizes the noise in the network and constructed a probabilistic Boolean net- 
work called Boltzmann machine. The deterministic and stochastic Boolean networks 
have found wide range of applications in biological networks. 

In our model, the states of the nodes(proteins, DNAs or RNAs) in the network 
at the n-th step are represented by variables = (X^,X^, ■ ■ ■ ,X^) respectively, 
where is the number of nodes in the network. Each node i has only two values, 

= 1 and = — 1, representing the active state and the reset state of this node 
respectively. 

Deterministic Boolean network model: 

The deterministic model consider in the present paper is a deformation of model 
Al in supporting information. Let us suppose A" is a fixed integer, 5" = {1, 2, • • • , A^}. 
We take the state space as {—1, 1}"^. Denote the state of the n-th step as X^ — 
{X^, Xl, ■ ■ ■ , X^), then the dynamic is as follows: 

If A:„7^ (-1,-1,-- - ,-l),then 
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1 x>0; 

where the function sign(x) = •{ x — 0; Hi — Ylf=i '^ij^l ~ is the input to 

-1 X < 0, 

the i-th node and Ui {1 <i < N), given a priori, are called the threshold of the i-th 
unit. 

And if Xn = (—1, —1, ■ ■ ■ , —1), then Xn+i = Xn when the parameter 7 = rep- 
resenting that the system is under normal environment, and Xn+i = (1,— 1,— I,-- - 1) 
when the parameter 7 = 1 representing there are some perturbation (signal) of the 
system (e.g. DNA damage) emerge. 

Stochastic Boolean Network model: 

The stochastic model consider in the present paper is a deformation of model 
C in supporting information, which is similar to the model in (Ge et al. 2007). 
Consider a Markov chain {X„ = (X^, X^, • • • , X^),n = 0, 1, 2, • • • } on the state 
space { — 1, 1}"^, with transition probability given as follows: 

N 

P(X„+l|XO=n^^(^n+ll^n), (2) 
i=l 

where 



in which Hi — Y^^=i '^ij^l ~ is the input to the i-th node, if Ylf=i '^ij^l Ui 
and Xn 7^ (-1, -1, • • • , -1) or i > 2, 



yi yi . 

1 — ^n' 



l+e-« ' n+1 

e~" -y-i _ -1 yi 

l+e-« ' n+1 -^ni 



if Ef=i TijXl = Ui and X^ ^ (-1, -1, • • • , -1) or i > 2, and 



1-7, K+i = K. 



if Xn = (— 1, — 1, • • • ,—1), where the parameter < 7 < 1 still represents the 
stochastic perturbations of the system from extracellular environment. 

In the above equation, a{> 0), j3{> 0), Tij and Ui (1 <i,j< N) are parameters 
of the model. The positive temperature-like parameter P represents noise in the 
system from the perspective of statistical physics (Amit 1989, Zhang et al. 2006). 
Noticeably, the actual noises within a cell might not be constant everywhere, but here 
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we use a system-wide noise measure for simplicity. To characterize the stochasticity 
when the input Hi to a node is zero, we have to introduce another parameter a. 
This parameter controls the likelihood for a protein to maintain its state when there 
is no input to it. 

The previous parameters (3 and a represent the intracellular noise due to ther- 
modynamic fluctuations, and on the other hand, we need another parameter 7 to 
characterize the extracellular signal strength with appropriate stochasticity. Since 
the signal is not purely disordered, we can not express the probability Pi(X^^_]^|X„) 
when Xn = (— 1,— I,-- - ,—1) as the exponential form analogous to the Boltz- 
mann distribution, hence we directly regard 7 as the probability transiting from 
,-1) to (1,-1,- ■■ ,-1). 

In our models below, Tij = 1 for a arrow from protein j to protein i, and 
Tij = — 1 for a horizontal bar instead of arrowhead from protein j to protein ^(Fig. 
[T],[2]). And it is indispensable to point out that self-connections haven't been taken 
into consideration in our models for simplicity. 

Similar to Proposition II. 6 in supporting information, with the same initial dis- 
tribution, when 7 — > 0, and a, (3 00, then the stochastic Boolean network model 
converges to the corresponding deterministic Boolean network model with parame- 
ter 7 = 0; and when 7 — > 1, and a, /5 — » 00, then the stochastic Boolean network 
model converges to the corresponding deterministic Boolean network model with 
parameter 7 = 1. 

Ill Steady states and synchronized dynamics of 
the p53 pathways 

From biochemical perspective, the microscopic variables for a cellular regulatory 
network are the concentrations, or numbers, of various mRNAs, regulatory pro- 
teins, and cofactors. If all the biochemistry were known, then the dynamics of such 
a network would be represented by a chemical master equation (McQuarrie 1967, 
Gillespie 1977). Unfortunately, much of the required information is not available, 
nor such a "fully- detailed" model will always be useful. Phenomenologically the 
concentrations of key players of a biochemical regulatory network can often be re- 
duced to two or three states, such as resting state, activated state, inactivated state. 
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etc. (Fall et al. 2002). The interactions between these states are usually determined 
from experimental data. 

The present study will build simple deterministic and stochastic Boolean network 
models for several negative feedback loops of the p53 pathways, and more impor- 
tant, provides a sound mathematical explanation of the synchronized dynamics and 
stochastic limit cycles in the stochastic Boolean network model after being perturbed 
by stress signals, in terms of the theory of nonequilibrium circulations (Jiang et al. 
2004). This makes the description of pathway more definite and penetrating. 

III.l The core regulation 

Negative feedback loops, composed of one transcription arm and one protein-interaction 
arm, are a common network motif across organisms. p53, the tumor suppressor, 
transcriptionally activates Mdm2, which in turn targets p53 for degradation. Al- 
though it is believed that the existence of negative feedback loop in biological systems 
can always gives rise to oscillations, yet it is not sufficient for the appearing of a 
limit cycle in the deterministic ODE model from the mathematical point of view. 

A bit out of expectation is, if we use the deterministic Boolean network even 
to the simplest case(FigIl]), the existence of negative feedback already gives rise to 
a limit cycle corresponding to oscillations in a biological system, when there exists 
the outside signal (i.e. 7 = 1). Fig 1. 

Model: From Section II(For FigH]). 

The first node X^: p53; and the second node X^: Mdm2. 

= 2 is the number of nodes in the model. For simplicity, we set thresholds as 
Ui = U2 = 0. And the interacting matrix 



T 



-1 

1 



The deterministic model with 7 = has a global attractor (—1,-1), which 
corresponds to the fact that the main function of the negative feedback between p53 
and Mdm2 is to keep p53 at a low steady state level in normal cells. 

On the other hand when 7 = 1, the deterministic model has a unique limit cycle 
consist of 4 state, which is described by Table [H It corresponds to the fact that the 
stress signal (e.g. DNA damage) will activate the protein p53(i.e. time-2 point in 



7 



Table [T]) and induce a transition to oscillations of p53 level after being perturbed 
from the outside environment. 

Time p53 Mdm2 



11-1 

2 1 1 

3 -1 1 

4 -1 -1 



Table 1: Cycle evolution of the protein states in the negative feedback loop of p53 
and Mdm2 after being perturbed. 

Moreover, as long as the stress signal strength 7 is sufficiently low (0 < 7 << 1) 
or sufficiently high (0<1 — 7<<1), the stochastic Boolean network model would 
preserve the dynamics of the corresponding deterministic model at a certain low 
level of noise. As it illustrates the same phenomenon as the Cyclin G/Mdm-2 loop 
given below, so we won't show the data here. 

Similar model can also be built and analyzed exactly by the same reasoning, to 
other ubiquitin ligases that promote p53 ubiquitination and subsequent proteasomal 
degradation (Fig. 10 in Harris et al. 2005). We only use the simplest example as a 
start, and pass directly to a really delicate interesting case. 



III.2 Cyclin G/Mdm2 loop 

Model: From Section II. 

The first node X^: p53; The second node X^: cyclin G; The third node X^: 
PP2A cyclin G; The last node X^: Mdm-2. 

N=4 is the number of nodes in the model. For simplicity, we set thresholds as 
Ui = U2 = U3 = U4 = 0. And from Fig|2l the interacting matrix 



T 



-1 
10 
10 
10 10 



The deterministic model with 7 = has a global attractor (— 1, — 1, — 1, — 1), 
which corresponds to the fact that p53 is kept at a low steady state level in normal 



Fig 2. 
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cells. 



On the other hand, the deterministic model when 7 = 1 has a unique limit cycle 
consist of 8 state, which is described by Table [2l It corresponds to the fact that the 
stress signal (e.g. DNA damage) will activate the protein p53(i.e. time-2 point in 
Table [2]) and induce a transition to oscillations of p53 level after being perturbed 
from the outside environment, which roughly corresponds to the p53 pathway or 
biological trajectory described in (Harris et al. 2005) "One of the most active of the 
p53-responsive genes is the cyclin G gene. It is rapidly transcribed to high levels 
after p53 activation in a wide variety of cell types." (i.e. Time-3 point in Table [2]) 
"The cyclin G protein makes a complex with the PP2A phosphatase, which removes 
a phosphate residue from Mdm-2, which is added to the Mdm-2 protein by a cdk 
kinase" (i.e. Time-4,5 points in Table [2]). Our language in the present paper is more 
definite and penetrating than the quotations. 



Table 2: Cycle evolution of the protein states in the Cyclin G/Mdm-2 loop after 
being perturbed. 

Numerical Simulation of the stochastic Boolean network model 

The numerical results about the cyclic motion of this stochastic model are quite 
similar to the stochastic Boolean network model of the cell-cycle, which we have 
recently discussed (Ge et al. 2007). The main conclusion is that, given the struc- 
ture of the negative feedback present in Fig. [2], the stochastic model still exhibits 
well pronounced oscillations after being seriously perturbed, which is excellently 
characterized by the circulation theory introduced in the supporting information. 

Numerical computations for the current model(Fig|2]), are carried out with the 



Time p53 Cyclin G PP2A Cyclin G Mdm-2 
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famous Gillespie's method (Gillespie 1977) of the stochastic Boolean network model 
using MATLAB, and the results are given in the following figures. The network 
with 4 binary nodes has a total of 2^ = 16 number of states. Here, we can present 
the dynamics of the network in terms of the integer states 0, 1, 2, ■ ■ ■ , 2^ — 1 = 15 
on a line. This 1-d system is reversible if and only if the 4-d system is reversible. 

Fig. [3] are the basic behavior of a random trajectory. The upper panel shows 
that there arises the phenomenon of local rapid synchronization like that observed 
in (Hopfield 1995) during a very short time period after the value of 7 transits 
from 0.01 to 0.99 at time n = 50, when jS and a are sufficiently large. The lower 
panel is a random trajectory over a longer time. Little deviation is shown from the 
deterministic trajectory in Table [2] after being perturbed at time n = 50, which 
implies that the stochastic model still leads to well pronounced oscillations when 
the perturbation from the environment is sufficiently high. Fig 3. 

Fig. m shows the stationary distribution of the state (—1,-1,— 1,-1) in the 
stochastic model which increasingly approaches 1 when f3 tends to infinity. This ex- 
cellently corresponds well to the dynamics of the deterministic model when 7 = 0. 
At large (3 (low temperature or small noise level), the low level state (—1, —1, —1, —1) 
is the most probable state of the system. So analogous to the concept of the deter- 
ministic model, this state (—1, —1, —1, —1) can be regarded as the global attractor 
of the stochastic model. Moreover, one observes a phase-transition like behavior 
of the stationary distribution of the state (—1,-1,— 1,-1) while varying the noise 
level f3 (similar behavior has been seen in (Qu et al. 2002)). Fig 4. 

Then we turn to investigate the synchronized dynamics when the signal pa- 
rameter 7 is sufficiently high(i.e. < 1— 7 << 1). The keys to understand 
synchronization behavior in stochastic Boolean network models are (i) establishing 
a correspondence between a stochastic dynamics and its deterministic counterpart; 
and (ii) identifying the cyclic motion in the stochastic models. 

As there is a growing awareness and interest in studying the effects of noise in 
biological networks, it becomes more and more important to quantitatively char- 
acterize the synchronized dynamics mathematically in stochastic models, because 
the concepts of limit cycle and fixed phase difference no longer holds in this case. 
Instead, physicists and biologists always have to characterize synchronized dynam- 
ics by the distinct peak of some spectrum or just only by observing the stochastic 
trajectories, which however may cause ambiguities in the conclusion. Therefore, 
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a logical generalization of limit cycle to the stochastic model is well worth to be 
developed. 

In case of the stochastic models of biological networks, there does exist a rather 
complete mathematical theory for the cyclic motion of the corresponding Markov 
chains, which has been developed for more than twenty years (Jiang et al. 2004, 
Kalpazidou 1995). One of the most important concepts in this mathematical theory 
of NESS is the circulation, which corresponds to the cycle kinetics in open chemical 
systems (Hill 1989, Qian 2005). The details of mathematical theory is supplied in 
the supporting information. 

To further characterize the synchronized dynamics, we give Fig. [5] which shows 
the Fourier power spectrum of the stochastic trajectory with different values of 
j3 respectively. Using MATLAB, the discrete Fourier transform for time series 
{xi, X2, - ■ ■ , Xn} is defined as 

m — 1 



^. ^-»(27r/n)(m-l)(fc-l) 



k=l 



-27r, l<m<n \ . (3) 



n 

Therefore, by the Herglotz theorem (Qian et al. 1997 p. 331), the power spectrum 
of discrete trajectory has a symmetry Um = yn+2-m- For different sets of parameters, 
we found all the calculations give the same outstanding main peak in the Fig. [51 It 
is important to mention here that by ergodicity, different trajectories give the same 
power spectrum for any (3 that is sufficiently large. Fig 5. 

The single dominant peak in Fig. [5] implies there exists a global synchronization 
and a globally attractive phenomenon. Note that by representating our maps, one- 
to-one, from the binary nodes to the integers — 15, the synchronized behavior 
is preserved. It is possible that the map will cause some distortion in the power 
spectrum. 

To further illustrate the synchronized behavior. Fig. [6] shows the power spectra 
of all the 4 individual nodes in the network. While subtle details are different, all ex- 
hibit the dominant peak, similar to that of the overall dynamics. This demonstrates 
further that synchronized dynamics is presented in the network. Fig 6. 

Fig. [7] plots the magnitude and the period of the dominant peak of the power 
spectrum in Fig. [5] respectively, as functions of the noise strength, i.e., the parameter 
j3. It shows, as we have predicted, that the period converges to 8 which corresponds 
perfectly to the number of states in the main cycle of Table [2] when (3 is large. We 
also put error bars on the upper panel of Fig. [7] with various values of 13. Fig 7. 
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Finally, Fig. [8] shows how the net circulation of the dominant, main cycle varies 
with applying the determinant presentation of circulations according to Theorem 
II. 4 in the supporting information. Note that circulation is just the time-averaged 
number of appearance of certain cycle along the stochastic trajectory of our model, 
and net circulation is just the difference between the circulations along the positive 
direction and negative direction of a specific cycle. It is clearly seen that the net 
circulation of the main cycle increases monotonically to | that is just the reciprocal 
of the number of states in the main cycle, which implies the appearing of more 
and more distinct synchronization and global attractive behavior with increasing 
(3. The direction of the net circulation does not change. It is also found that the 
circulation of negative direction along the main cycle is always quite low compared 
to the circulation of positive direction. Fig 8. 

The net circulations of all the cycles are very small when j3 and a are near zero, 
since the system is close to equilibrium (reversible) state when a = (3 = Q according 
to lemma II. 9 in supporting information. 

For large /5, the net circulations of all the cycles except the main cycle are also 
very small by numerical simulation using the determinant expression in Theorem 
II. 4 of supporting information. All the circulations of non-main cycles actually de- 
crease with increasing (3 when (3 is large. Examples are shown in Fig. [91 This large 
separation in the magnitudes of the weights gives rise to the stochastic synchro- 
nization, and this stochastic limit cycle can be defined as a "stable" one, whose 
attractive domain is global. Fig 9. 

As in (Zhang et al. 2006, Ge et al. 2007), we also notice that there exists an 
inflection point in the curve in Fig. [HI This implies a cooperative transition of 
the net circulation of the main cycle while varying the noise level /3, which equiv- 
alently means some sort of "phase transition" around (3 = 2 from chaotic fiuctu- 
ations(period=infinity) at small (3 to periodic fluctuations (period=8) at large [3. 
Certainly, the implication of this observation remains to be further elucidated. 

IV Discussion 

From detailed models to simplified Boolean network approach 

Deterministic nonlinear mathematical models, based on the Law of Mass Action, 
have been traditionally used for biochemical reaction networks. Furthermore, noises 
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are unavoidable in small biochemical reaction systems such as those inside a single 
cell. Stochastic models with chemical master equations (CME) (McQuarrie 1967, 
Van Kampen 1981) should be developed, which has already provided important 
insights and quantitative characterizations in some cases of the biochemical system 
(Fox et al. 1994, Fox 1997, Zhou ct al. 2005). Ref. (Wilkinson 2006) is a good 
introduction to the stochastic modelling in biology. 

In many cellular biochemical modelling, it is impossible to build a detailed, 
molecular number-based CME model due to a lack of quantitative experimental 
data. Thus, one often seeks a simplified network model based on simple binary 
states of the signaling molecules. This leads to the Boolean, Markov network model 
of the p53 pathways we studied in the present paper. 

Modelling the oscillatory dynamics of p53 pathways after being per- 
turbed 

The realization that p53 is a common denominator in human cancer has stim- 
ulated an avalanche of research since 1989. The p53 gene can integrate numerous 
signals that control cell life and death, and damped oscillations for p53 and Mdm2 
has been observed and modeled (Lev Bar-Or et al. 2000, Ma et al. 2003). The 
Boolean network models built in the present paper, which omit some parameters 
to represent the degradation of the p53 protein in cells, only try to explain the 
mechanism of the oscillatory behavior of the p53 dynamics after being perturbed by 
stress signals rather than exhibiting the damped oscillations. The signal parameter 
7 plays a very important role in the model similar to the cell cycle model in (Ge et 
al. 2007), which can induce the level of p53 from a low steady state to oscillated 
behaviors. 

On the other hand, ideally, one should try to combine these Boolean network 
models of negative feedback loops together to construct a clear and integrated picture 
of the p53 pathways, maybe especially including several positive feedback loops 
(Harris et al. 2005). But unfortunately we haven't developed a reasonable way to 
do so. 

Synchronization and circulation in stochastic Boolean network 

Synchronization is an important characteristics of many biological networks 
(Strogatz 2003, Winfree 2000) whose dynamics has been modelled traditionally by 
deterministic, coupled nonlinear ordinary differential equations in terms of regu- 
latory mechanisms and kinetic parameters (Murray 2002, Fall et al. 2002). Two 
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important classes of biological networks which have attracted wide attentions in 
recent years are neural networks (Scott 2002) and cellular biochemical networks 
(Goldbeter 1996). 

As we know, the occurrence of a deterministic limit cycle in an ODE model is 
the hallmark of a synchronization phenomenon, while such a definite concept no 
longer holds in a stochastic system. It is observed that in our present model of p53 
pathways, the trajectory concentrates around a main cycle, which we call stochastic 
limit cycle and is very similar to in many aspects the limit cycle of a deterministic 
model. Furthermore, the present work shows that the circulation in an irreversible 
Markov chain (Jiang et al. 2004) is a reasonable generalization of the concept of 
deterministic limit cycles, and provides a sound mathematical explanation of the 
synchronization in the stochastic Boolean network model. 

Stability and robustness of the p53 genetic networks 
Biological functions in living cells are controlled by protein interaction and ge- 
netic networks, and have to be robust to function in complex (and noisy) environ- 
ments. More robustness would also mean being more evolvable, and thus more likely 
to survive. 

More precisely, these molecular networks should be dynamically stable against 
various fluctuations which are inevitable in the living world. Therefore, the stability 
and robustness of Boolean network model can not be determined if the noise hasn't 
been introduced into the model, since it is not reasonable to simply apply the Eu- 
clidean topology in mathematics to such a real biological discrete model not to say 
the magnitudes of the different attractive basins of fixed states and limit cycles as 
in the deterministic models. For instance, the recent analysis of stochastic cell-cycle 
Boolean networks (Ge et al. 2007, Zhang et al. 2006) just claimed the stabihty and 
robustness of the previous deterministic model (Li et al. 2004), in which Li, et.al. 
announced but didn't really investigate the robustness of this model against pertur- 
bation. Hence, the numerical and theoretical studies in the present paper in some 
sense excellently exhibit the stability and robustness of the p53 genetic networks. 
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Figure 1: Negative feedback loop of p53 and Mdm2. Arrows denote stimulatory 
interactions, whereas horizontal bars instead of arrowheads indicate inhibitory in- 
fluences. 
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Figure 2: Cyclin G/Mdm2 loop, redrawn from (Fig. 8 in Harris et al. 2005). 



Local rapid synchronization 




Figure 3: Stochastic trajectory and synchronization of the Cyclin G/Mdm-2 loop. 
Simulations are carried out with the parameters a = 5, f3 = 5 and the values of 7 
transits from 0.001 to 0.999 at time n = 50. 
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Figure 4: Stationary distribution of the state 
model. Simulations are carried out with the parameters a = 5 and 7 = 0.001. 
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Figure 5: Power spectrum of the overall trajectory in the Cyclin G/Mdm-2 loop 
with the parameter a = 5 and 7 = 0.999 fixed, and with different (3: blue: f3=2A] 
green: /3=4.8; red: f3 = Q. The discrete Fourier transform causes an alias; hence the 
spectrum is symmetric with respect to tt on the [0, 27r] interval. 
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Figure 6: Power spectra of individual nodes in the Cyclin G/Mdm-2 loop show a 
synchronization among all the nodes with the parameter o; = 5 and 7 = 0.999 fixed, 
and different j3\ blue: /5=2.4; green: /3=4.8; red: /3=6. 
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Figure 7: Magnitude and period of the dominant power-spectral peak as functions of 
P in the Cychn G/Mdm-2 loop, with a = 5 and 7 = 0.999. In the upper panel, the 
magnitude of the dominant power-spectral peak is averaged over 20 simulations. The 
solid curve is the mean, and the dotted curves are the mean ± standard deviation. 
The period approaches to 6 (the horizonal line) with increasing (3. 
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Figure 8: Net circulation of the main cycle in the Cyclin G/Mdm-2 loop as a function 
of the noise strength, ^, with the parameter a = 5 and 7 = 0.999 fixed. 
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Figure 9: Net circulation of several other cycles in the Cyclin G/Mdm-2 loop as a 
function of the noise strength, /3, with the parameter a = 5 and 7 = 0.999 fixed. 
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Deterministic and Stochastic Boolean Networks 

This section is mainly restated from (Ge et al. 2007). 



> 
(N 

vn 

<N ; We first give a brief account of the deterministic Hopfield networks and its 

Q>^ ■ stochastic incarnation, the stochastic Boolean networks (sometimes called Boltz- 

o ■ 

mann machines). They can be mainly categorized into several classes as follows 

X 



^ ; (Amit 1989). 



Let us suppose is a fixed integer, S = {1, 2, ■ ■ ■ , A^}. We take the state space 
as {-1,1}^. 

Model A: deterministic 

Al (discrete time, synchronous, McCulloch-Pitts): Denote the state of the n-th 



*School of Mathematical Sciences, Peking University, Beijing 100871, P.R.C.; email: ed- 
mund_ge@tom.com 

^School of Mathematical Sciences, Peking University, Beijing 100871, P.R.C. 



1 



step as Xn = (X^, X^, ■ ■ ■ , X^), then the dynamic is 

(N \ N 

Y,T.,xi-uA, i^Y,T^,xi^u,■ (1) 
i=i / j=i 

and if X]jLi^ij"^n ~ then X*^^ randomly choose 1 or —1 with probabihty ^ 
respectively. Ui {1 < i < N), given a priori, are called the threshold of the ith unit. 

A2 (continuous time, synchronous): Every state has an exponentially distributed 
stochastic waiting-time, with mean waiting-time A~^, then chooses the next state 
by the same rule of model Al. 

A3 (discrete time, asynchronous, Hopfield): The neurons are updated one by 
one, in some prescribed sequence, or in a random order. If the previous state cr 
satisfies J2f=i'^ij^j > then cxj changes to be 1, otherwise changes to be —1. 

A4 (continuous time, asynchronous): Every state has an exponentially dis- 
tributed stochastic waiting-time, with rate constant A, then chooses the next state 
by the same rule of model ^43. 

Remark I.l Note that by deterministic, we mean the transition from one state to 
the next is deterministic. But the systems with continuous-time will still behave 
stochastically due to the Poisson nature in the transition time. 

Model B: stochastic Boolean networks 

Bl (discrete time, synchronous): Consider the Markov chain 

{X„ = {Xl X^, . . . , X„^), n = 0, 1, 2, ■ ■ ■ } (2) 



"'^The function sign(x) = < 



1 x>Q] 
a; = 0; 
-1 a; < 0. 



on state space {—1, 1}'^, with transition probability given as follows: for each pair 
of states (J,ri E { — 1, 1}'^, the probability transiting from a to r] 

where /?(> 0), and Ui {1 < i, j < N) are parameters of the model. 

B2 (continuous time, synchronous): Consider the continuous-time Markov chain 
{^t '■ t ^ 0} on state space {—1, 1}'^. Every state waits an exponential time with 
meantime until choosing the next state by the rule of model Bl. So the transition 
density matrix is 

qan^>^Pav: V (7, 77 G {-1, 1}^. (4) 

B3 (discrete time, asynchronous): Denote cr* to be the new state which changes 
the sign of the ith coordinate of a. Consider the Markov chain {X^ = {X^, X^, • • • , ) , n 
0, 1, 2, • • • } on state space {—1, 1}'^. The neurons are updated one by one, in some 
prescribed sequence, or in a random order. Then choose the next state according to 
the probability: 

^ exp(/3a4Ej^i T^j-gj - Uj)) ^ ^ |_i ns 

^""^ ~ exp(/3(E7=i 71,<7, - [/,)) + exp(-/?(E7=i T,,<7, - C/.)) ' 

where P > 0; and = 0, if 77 7^ cr* for each i. 

B4 (continuous time, asynchronous): Every state waits an exponential time 

with meantime until choosing the next state by the rule of model B3. So the 

transition density matrix is 

Qaa^ = ^Paa-, VcT G {-1, l}^. (5) 

The third class given below is a variant of the model Bl. It is included here since 
it is the model used for the probabihstic Boolean network of cell-cycle regulation in 



(Geet al. 2007). 
Model C: deformation 

Fix a > 0. Consider a new Markov chain {X„ = (X^^X^,-- - ,X^),n — 
0, 1, 2, • • • } on the state space {—1, 1}'^ taking the model Bl as defined initially, 
with transition probability given as follows: 

N 

P(X„+i|X„)=n^.(^;+ll^n), (6) 

1=1 

where we define 



Pi{Xl^_^^\Xn) - 



exp(/3(Ef=i rij^^-;^i))+exp(-/3(Ef=i Ti^Xi-Ui)) 



l+e 



if Y^U^an-U.^ndXi^.^Xl 



l+e-« 



if Y.U^i,Xi = U, and = l-Xl 



Remcirk 1.2 The model C differs from Bl when Yl!j=\ '^ij-^l — Ui = 0, and the 
latter is also a special case of the former when a — 0. 



II Mathematical Circulation Theory of Network 
Nonequilibrium Steady States 

This section is also mainly restated from (Ge et al. 2007), ignoring those detailed 
proofs. 

There are many different approaches to the theory of nonequilibrium statistical 
mechanics in the past (Nicolis et al. 1977, Keizer 1987), mathematical theories of 
which have emerged in the last two decades, and Jiang et.al. have summarized 
their results of this theory in a recent monograph (Jiang et al. 2004). The most 
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important concepts in the theory are (i) reversibihty of a stationary process that 
corresponds to thermodynamic equihbrium, and (ii) the circulation in a stationary 
process which corresponds to NESS. A key result of the theory is the circulation 
decomposition of NESS. 

II. 1 Circulation theory of nonequilibrium steady states 

Hill (Hill 1989) constructed a theoretical framework for discussions of vivid metabolic 
systems, such as active transport, muscle contractions, etc. The basic method of 
his framework is diagram calculation for the cycle flux on the metabolic cycles 
of those systems (Hill 1989). He successively found that the result from diagram 
calculation agrees with the data of the numbers of completing different cycles given 
by random test (Monte Carlo test), but did not yet prove that the former is just 
the circulation rate in the sense of trajectory of a corresponding Markov chain. 
In Chapter 1 and 2 of (Jiang et al. 2004), Markov chains with discrete time and 
continuous time parameter are used as models of Hill's theory on circulation in 
biochemical systems. The circulation rate is defined in the sense of trajectories and 
the expressions of circulation rate are calculated which coincide with Hill's result 
obtained from diagrams. Hence the authors verify that Hill's cycle flux is equivalent 
to the circulation rate defined in the sense of trajectories. 

Below we only state the circulation theory of discrete-time Markov chains, and 
refer to Chapter 2 in (Jiang et al. 2004) for the quite similar circulation theory of 
continuous-time Markov chains. 

First of all, we state the main results, rewritten in terms of the cycle representa- 
tion of stationary homogeneous Markov chains (Kalpazidou 1995, Theorem 1.3.1), 
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which is analogous to the Kirchoff 's current law and circuit theory in the networks 
of master equation systems (Schnakenberg 1976). 

Theorem II. 1 Given a finite oriented graph G = {V,E) and the weight on every 
edge {we > : e G E}. If the weight satisfies the balance equation (input= output) 

J] We = 5^ w,, Vz G V, (7) 

then there exist a positive function defined on the oriented cycles {wc '■ c = (ei, 62, • ■ ■ , Cfc), 
Z^}, such that 

We = ^w^, ye EE. (8) 

For a finite stationary Markov chain X with finite state space 5* = {1, 2, ■ ■ ■ , A^}, 
transition matrix P = {pij : i,jE S} and invariant distribution vf = {tti, 7T2,--- , tt^}, 
one can take its state space to be the group of vertexes, and E = {e : e~ = 
i,e+ = j,pij > 0} with weight {we = f^iPij}- Then notice that ([7]) is satisfied, since 

"P^i ~ ^ ea:Ch. i and 

i j 

we can conclude that 

Corollary II. 2 (cycle decomposition) For an arbitrary finite stationary Markov 
chain, there exists a positive function defined on the group of oriented circuits 
{wc : c = (ii, i2, ■ ' ' k E Z~^} such that 

T^iPij = ^WcJciiJ), yij e V, (9) 

c 

where Jc{i,j) is defined to be 1 if the cycle c includes the path i j , otherwise 0. 



Definition II. 3 Wc is called the circulation along cycle c. 
For any ij e S,ij^ j, 

TCiPij - TTjPji = ^(wc - WcJJciiJ), (10) 

c 

where c_ denotes the reversed cycle of c. Equation ffTOl) is called the circulation 
decomposition of the stationary Markov chain X, and the difference between 
the circulations of positive direction and negative direction along the cycle c (i.e. 
Wc — Wc_) is called the net circulation of cycle c, which actually represents the 
true flux of this cycle. 

It can be proved that generally the circulation decomposition is not unique, i.e. 
it is possible to find another set of cycles C and weights on these cycles {wc\c G C} 
which also fit (fTOl). 

However, the most reasonable choice of circulation definition is the one defined 
in the sense of trajectories form the probabilistic point of view. Along almost every 
sample path, the Markov chain generates an infinite sequence of cycles, and if we 
discard every cycle when it is completed and at the meantime record it down, then 
we can count the number of times that a specific cycle c is formed by time t, which 
we denote by Wc,t{^)- 

The following theorem is recapitulated from Theorem 1.3.3 in (Jiang et al. 2004). 

Theorem II. 4 Let Cn{uj), n = 0, 1, 2, ■ ■ ■ , be the class of all cycles occurring until 
n along the sample path {Xi{uj)}. Then the sequence {Cn{uj),Wc,n{^^)/n) of sam- 
ple weighted cycles associated with the chain X converges almost surely to a class 
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(Coo, Wc), that is, 

Coo = lim Cn{uj), a.e. (11) 

n— >+oo 

Wc = lim ^'^'"^ ^ , a.e. (12) 

Furthermore, for any directed cycle c = {ii,i2, ■ ■ ■ ,is) ^ Coo, the weight Wc is given 
by 

'^c — Piii2Pi2i3 ' ' ' Pis-lisPisil v-^ ^^^^ ■^r\ ' ^ ' 

where D = {dij} = I — P = {6ij — Pij} and D{H) denotes the determinant of D 

0, t^j; 



with rows and columns indexed in the index set H. The function 6ij 
is the well known Kronecker delta function. 



1, i = J, 



It is important to emphasize that the circulations defined in the above theorem 
also satisfy the circulation decomposition relation flTUl) . 

The above theorem not only rigorously confirms the Hill's theory, but also gives 
a prior substitute method of the widely used diagrammatic method. The complexity 
of directed diagrams and cycles increases rapidly with the number of states in the 
model, while the determinant interpretation is much more systematic and easy to 
be applied using the mathematics softwares. But it will still cost excessive time to 
compute these determinants if there are hundreds of states in the model. 

Luckily, a Monte Carlo method using the so-called derived chain method 
(Section 1.2 in Jiang et al. 2004) to compute the cycle fiuxes has already been 
developed according to the above theorem, the main idea of which is just simply to 
discard every cycle when it is completed and at the meantime record it down so as 
to count the number of times that a specific cycle c is formed by time t{i.e. Wc^ti^))- 
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Then when the time t is long enough, one gets the approximated circulation of the 
cycle c (i.e. Wc ~ 

Recently, the Boolean yeast cell-cycle network model discussed in (Ge et al. 
2007) has 2048 states and it is found that the Monte Carlo method is much more 
efficient than the method of determinant interpretation, because the latter is even 
impossible to be applied to such a large model upon a normal computer. 

The relationship between circulation and NESS is as follows, which is recapitu- 
lated from Theorem 1.4.8 in (Jiang et al. 2004). 

Theorem II. 5 Suppose that X is an irreducible and positive-recurrent stationary 

Markov chain with the countable state space S, the transition matrix P = ipij)ijes 

and the invariant probability distribution U — {i^iji^s, o-nd let {wc '■ c e C^} be the 

circulation distribution of X, then the following statements are equivalent: 

(i) The Markov chain X is reversible. 

(a) The Markov chain X is in detailed balance, that is, 

(Hi) The transition probability of X satisfies the Kolmogorov cyclic condition; 

Piii2Pi2i3 ' ' ' Pis-iisPish ~ PhisPisia-i ' ' ' PiahPhh t 

for any directed cycle c — (ii, • • • , i^). 

(iv) The components of the circulation distribution of X satisfy the symmetry con- 
dition: 

Wc = Wc_,\/c e Coo- 
Consequently, when the system is in a nonequilibrium steady state, there exists 
at least one cycle, containing at least three states, round which the circulation rates 
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of one direction and its opposite direction are asymmetric (unequal), so as to cause 
a net circulation of the cycle. In theoretic analysis, if there is a large separation in 
the magnitude of the circulation, between few dominant, main cycles and the rest, it 
gives rise to the stochastic synchronization phenomenon and helps to distinguish the 
most important main biological pathways, which can be observed in experiments. 

II. 2 Applied to deterministic and stochastic Boolean net- 
works 

The keys to understand synchronization behavior in stochastic models are {i) es- 
tablishing a correspondence between a stochastic dynamics and its deterministic 
counterpart; and (ii) identifying the cyclic motion in the stochastic models. 

In the framework of the stochastic theory, deterministic models are simply the 
limits of stochastic processes with vanishing noise. This is best illustrated in the 
following proposition. 

Proposition II. 6 With the same initial distribution, when (3 oo, the model Bk 
converges to the model Ak in distribution, for k = 1,2, 3, 4. 

We shall further discuss the necessary condition for stochastic Boolean networks 
to have cyclic motion. In the theory of neural networks, one of Hopfield's key results 
(Hopfield 1982,1984) is that for symmetric matrix Tij, the network has an energy 
function. In the theory of Markov processes, having a potential function (Kolmogorov 
cyclic condition in Theorem III.Sp is a sufficient and necessary condition for a re- 
versible process. A connection is established in the following theorem. 



10 



Theorem II. 7 For k — 1,2,3,4, the Markov chain {X^} in the model is re- 
versible if and only if T^j = Tji, Wi, j = 1,2, ■ ■ ■ , N , i.e. the matrix T is symmetric. 

Prom the above two conclusions, it is obvious that 

CorollEiry II. 8 When T^j — Tji,\fi,j — 1,2, - ■ ■ ,N, there doesn't exist any limit 
cycle consist of more than two states in model Ak, k = 1, 2, 3, 4. 

Lemma II. 9 For k — 1, 2, 3, 4, the Markov chain {Xn} in the model is reversible 
when (3 is zero, and the Markov chain {X^} in the model C is reversible when (3 
and a are both zero. 
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